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We  have  developedja  self-consistent  quantitative  model  to  compute  electric 
fields,  currents  andCthe  resulting  plasma  flow  in  the  inner-magnetosphere/ 
ionosphere  system  (L  010) ;  parallel  electric  fields  and  ionospheric  neutral 
winds  are  not  included.  The  model  was  tested  for  a  substorm-type  event  that 
occurred  on  September  19,  1976.  Satellite  data  (primarily  from  the  Air  Force 
S3-2  satellite)  were  used  extensively  both  for  boundary  conditions  and  for 
comparisons  with  model  predictions.  Other  data  were  also  used  as  input  for 
our  time-dependent  magnetic  field  and  conductivity  models. 

The  S3-2  data  for  the  event  show  some  novel  features,  independent  of  the 
simulation.  Dawn-dusk  electric  fields  show  a  general  correlation  with  the 
east-west  magnetic  field  perturbations.  Unexpectedly,  two  of  the  passes  dis- 
play  substantial  regions  of  sunward  plasma  flow  poleward  of  the  main  part  of  ei 
region-1  Birkeland  currents.  I 

The  cross-polar-cap  potential  drops  computed  from  the  data  represent  the 
first  effort  at  satellite  monitoring  of  this  important  parameter  during  various 
phases  of  a  substorm,  and  show  an  important  enhancement  during  the  substorm. 

Numerical  results  from  these  first-try  simulations  were  consistent  with 
most  of  the  established  features  of  convection  in  the  inner  magnetospherefc-suc 
as  generally  sunward  flow,  shielding  of  the  potential  electric  field  for  L  ^  b 
and  the  tendency  for  stronger  electric  fields  on  the  dusk  side  than  on  the  dawn 
side.  In  addition,  the  model  reproduces  some  typical  substorm  phenomena,  such  as 
energy-dependent  particle  injection  with  a  dawn-dusk  asymmetry  and  establishment 
of  a  partial  ring  current. 

This  paper  deals  with  model  logic,  methodology,  inputs  and  overview  of 
|  results;  the  succeeding  two  papers  give  detailed  analyses  and  comparisons  with 
I  data. 
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I.  INTRODUCTION 

The  present  three  papers  are  the  latest  reports  on  a  longstanding 
effort  at  a  quantitative  model  of  plasma  flow,  particle  distributions  and 
electric  currents  in  the  inner  magnetosphere  and  ionosphere,  and  at 
theoretically  modeling  the  time  dependent  physical  processes  that  take 
place. 

In  the  past  few  years  our  group  has  developed  a  program  for 
self-consistent  computer  simulation  of  the  magnetosphere- ionosphere 
system  [Wolf,  1970;  Jaggi  and  Wolf,  1973;  Wolf,  1974;  Hare!  and  Wolf, 
1976].  In  those  papers  the  emphasis  was  put  on  "steady  state 
simulation",  where  our  model  was  tested  for  "average"  magnetospheric 
conditions.  Over  this  period,  our  work  progressed  to  include  more  and 
more  physical  processes  while  the  mechanics  of  our  computer  simulation 
were  improved.  Although  some  interesting  and  meaningful  insights  into 
the  magnetosphere  were  obtained  from  such  quiet-time  simulations,  the 
magnetosphere  is  seldom,  if  ever,  close  to  a  steady  state,  and  the  most 
exciting  phenomena  occur  during  disturbed  times.  We  have  therefore 
proceeded  to  simulate  a  substorm-type  event  and  to  compare  our 
theoretical  results  with  satellite  data. 

The  concept  of  magnetospheric  convection  was  first  formulated  in  the 
early  1960's  [Axford  and  Hines,  1961;  Dungey,  1961;  and  Cole,  1961]. 

Some  sort  of  effective  friction  between  the  solar  wind  and  the 
magnetosphere  was  seen  to  cause  this  olasma  circulation,  resulting  in 
antisunward  flow  in  the  outer  magnetosphere,  and  sunward  flow  in  the 
inner  regions.  (See  reviews  by  Axford  [1969]  and  Stern  [1977]).  The 
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role  of  ionospheric  conductivity  and  Birkeland  currents  in  regulating 
this  sunward  convection  was  partly  understood  in  the  mid-sixties 
[Karlson,  1963;  Fejer,  1964;  Block,  1966].  However,  the  interactions  are 
many  and  complex,  and  many  of  the  subsequent  studies  modeled  only  a 
portion  of  the  system.  One  group  of  models  concentrated  primarily  on  the 
ionosphere  [Volland,  1975;  Heppner,  1977;  Yasuhara  and  Akasofu,  1977; 
Nisbet  et  al.,  1978;  Nopper  and  Carovi llano,  1978,  1979;  K amide  and 
Matsushita,  1979a  and  b;  Gizler  et  al . ,  1979].  Another  group  attempted 
modeling  the  ring  current  injection  using  assumed  semiempirical  electric 
fields  [e.g.,  Mcllwain,  1974;  Roederer  and  Hones,  1974;  Konradi  et  al., 
1976;  Cowley,  1976;  Kivelson,  1976;  Ejiri  et  al.,  1977,  1978;  Smith  et 
al. ,  1979].  Such  models,  though  useful,  were  generally  limited  in 
scope.  The  more  comprehensive  efforts  at  quantitative  modeling  of  the 
complete  magnetosphere- ionosphere  system  were  obviously  more  complex 
[Wolf,  1970;  Swift,  1971;  Vasyliunas,  1970,  1972;  Jaggi  and  Wolf,  1973; 
Mal'tsev,  1974;  Wolf,  1974;  Hare!  and  Wolf,  1976]. 

The  complexity  of  the  system,  its  nonlinearity  and  its  strong 
time-dependence  (especially  during  disturbed  times)  make  quantitative 
modeling  difficult.  We  have  developed  an  approach  of  coupling  the 
available  time-dependent  input  data  with  the  solution  of  the  differential 
equations  governing  the  large  scale  motions  of  plasmas  in  the 
magnetosphere- ionosphere  system.  While  some  of  the  data  are  used  as 
input,  other  data  are  used  for  testing  our  model  predictions.  We  have 
chosen  to  model  the  substorm-type  event  of  September  19,  1976  whose  onset 
occurred  at  roughly  1000  UT.  Given  input  data  (mostly  from  the  S3-2 
satellite)  our  program  computed,  self-consistently,  the  time  evolution  of 
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electric  fields,  plasma  flow  velocities,  electric  currents  and  ring 
current  injection.  Of  course,  Ine  fact  that  the  critical  spacecraft 
observations  cannot  be  made  continuously  in  universal  time,  at  all 
l-values  and  local  times  in  the  magnetosphere,  means  that  our  model  input 
data,  based  mainly  on  interpolations  of  available  observations,  contain 
far  less  spatial  and  temporal  structure  than  was  actually  present  on  19 
September  1976. 

It  should  be  emphasized  that  in  this  paper  as  well  as  in  the  two 
following  papers  we  have  not  "fudged*'  our  input  assumptions  to  improve 
agreement  between  theoretical  predictions  and  observations.  We  consider 
this  to  be  a  "true  first  run"  simulation.  Analysis  of  these  initial 
results  will  enable  us  to  improve  our  initial  and  boundary  conditions  to 
yield  better  agreement  with  data  in  future  simulations.  We  will  refer 
to  the  following  two  papers,  Harel  et  al.  [1980]  and  Spiro  et  al .  [1980] 
as  Papers  2  and  3,  respectively. 

While  the  present  first-run  simulations  are  an  imprecise 
representation  of  what  actually  happened  in  the  magnetosphere  and 
ionosphere  on  19  September  1976,  they  provide  many  important  insights 
into  the  behavior  of  the  inner  magnetosphere  during  substorms  and  other 
periods  of  strong  convection. 
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II.  FORMULATION  AND  LOGIC 

Our  model  concentrates  on  the  Inner  magnetosphere,  specifically  the 
region  where  magnetic  field  lines  are  certainly  closed.  We  feel  that  the 
complexity  of  the  outer  magnetoispheric  dynamics  does  not  lend  itself  to 
detailed  quantitative  modeling  yet.  Therefore,  we  apply  our  boundary 
conditions  just  inside  of  the  magnetopause  boundary  layer,  which  roughly 
maps  to  the  polar-cap  boundary,  and  attempt  a  self-consistent  solution 
equatorward  of  that  region  (Fig.  1). 

A.  Conservation  equations 

The  coupling  between  the  magnetosphere  and  the  ionosphere  requires 
simultaneous  solution  of  the  current-  and  particle-conservation  equations 
in  both  regions.  Drifts  of  plasma-sheet  and  ring-current  particles  are 
computed  under  the  assumption  that  the  pitch  angles  of  the  particles  are 
scattered  many  times  in  the  time  it  takes  them  to  drift  a  significant 
distance  in  the  magnetosphere  [Vasyliunas,  1968;  Kennel  1969],  Within 
this  isotropic-pitch-angle  approximation,  it  is  convenient  to  define  a 
parameter  A,  which  is  related  to  the  thermal  kinetic  energy,  E,  by 

E  =  A(/  ds/B)'2/3  (1) 

where  Zds/B  =  volume  of  a  magnetic  flux  tube  with  unit  magnetic  flux.  It 
is  shown  in  the  Appendix  that  this  "energy  invariant".  A,  and  also  the 
"nunber  invariant",  n{ X) .  which  is  the  nunber  of  particles  per  unit 
magnetic  flux  with  energy  invariant  A,  are  constant  along  the  drift  path 
of  a  particle. 

We  characterize  the  plasma-sheet  particle  population  in  terms  of  21 
species,  such  that  each  particle  of  species  k  has  charge  q^  and  energy 


invariant  Ten  species  are  ions,  ten  are  electrons  and  one  is  cold 
plasma;  electron  precipitation  is  taken  into  account  approximately  by 
considering  two  species  for  each  A-value,  corresponding  to  different 
amounts  of  loss,  as  discussed  in  section  V.  For  simplicity,  we  assume  a 
uniform,  steady  particle  source  in  the  tail  so  that  n^x)  is  equal  to  a 
constant  value  n^.  (Earthward  of  the  inner  edge  for  species  s,  is, 

of  course,  zero). 

The  drift  velocity  of  the  equatorial  crossing  point  of  a  particle  of 
species  k  is  given  by 

tk  ‘  Be'  *  *  <Veff  ?fnd  >  <2> 

where 

Veff  1  V  ♦  Vcor  ♦:  (Vqk>(/ds/B)’2/3  (3) 

Bg  =  equatorial  magnetic  field,  z  =  unit  vector  normal  to  the  equatorial 
plane  (northward);  Vg  =  gradient  operator  in  the  equatorial  plane,  V  = 
electrostatic  potential  in  a  reference  frame  that  rotates  with  the  earth, 
VC()r  *  corotation  potential  [equation  3  of  Jaggi  and  Wolf,  1973],  E-nd  * 
induction  electric  field. 

We  have  chosen  E^  such  that 

~ecp  (Je*^  =  ('ind  (*e*t)xU/Be2  (4> 

when  v_  (x  , t)  is  the  velocity  of  the  equatorial  crossing  point  of  a 
magnetic  field  line  whose  equatorial  crossing  point  is  xg  at  time  t  and 
whose  ionospheric  crossing  point  is  constant  in  time. 
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The  third  term  in  equation  3,  which  represents  gradient  and 
curvature  drift,  is  derived  in  the  Appendix.  We  use  equation  2  to  follow 
the  motion  of  the  inner  edge  of  each  particle  species  k  (see  section  V.C). 

The  gradient/curvature  drift  current,  mapped  along  field  lines  to 
the  equatorial  plane,  is  given  by  , 

4e(;e)  *  j[Xknk  *  x  Ve  C(/ds/B)"2/3]  (5) 

where  j  =  current  per  unit  length  in  the  equatorial  plane,  and  the  sum 
includes  the  species  k  that  are  present  at  xg.  Conservation  of  current 
in  the  magnetosphere  then  gives 


=  -V  *j 
e  -e 


-  *Jz  xV 


•f 


[(/ds/B)'2/3] 


(6) 


where  jj|e  =  Birkeland  current  per  unit  area  away  from  the  equatorial 
plane.  Magnetization  current  has  been  ignored  becuase  it  has  no 
divergence.  Note  that  (6)  implies  that  Birkeland  currents  are  generated 
only  near  the  inner  edge  of  the  plasma  sheet,  where  there  are  gradients 
in  the  n^'s. 

We  now  equate  the  current  out  of  the  equatorial  plane  to  the  total 
current  into  the  ionosphere  (both  hemispheres).  Assuming  that  the 
induction  electric  field  is  very  small  compared  to  the  potential  electric 
field  in  the  ionosphere,  neglecting  vertical  electric  field  in  the 
conducting  region  of  the  ionosphere,  and  neglecting  vn  x  B,  where  ,yn  * 
neutral-wind  velocity,  we  can  derive  a  simplified  form  of  the 
ionospheric-current  conservation  equation: 
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v  •[£•  (-vhv)]  = 

h  = 


(7) 


where  represents  a  horizontal  gradient  operator  in  the  ionosphere,  |  = 
height-integrated  ionospheric  conductivity  (both  hemispheres),  and  , 
the  current  per  unit  area  down  into  the  ionosphere,  is  given  by 


J|1  ■  <|B1r|/Be»|le  (8) 

where  B^p  «  radial  component  of  the  magnetic  field  at  the  ionosphere.  By 
using  the  same  V  in  (3)  and  (7)  we  implicitly  assume  E  •  B  =  0. 

B.  Logical  loop 

Table  1  summarizes  the  assumptions  of  our  model;  model  logic  is 
illustrated  in  Figure  2.  The  basic  logical  loop  (the  central  pentagon  of 
the  figure)  is  a  modification  of  a  diagram  given  by  Vasyliunas  [1970], 
Dashed  lines  mark  future  additions  to  our  program  which  are  not  included 
in  the  present  simulation.  Starting  at  a  given  time  T  with  a  given 
hot-particle  distribution  in  the  magnetosphere  (section  IV. D),  estimated 
from  average  plasma-sheet  particle  data  in  the  present  case  (upper  box), 
we  proceed  counter-clockwise.  We  first  solve  equation  6  for  the 
divergence  of  the  gradient/curvature  drift  current  in  the  magnetosphere. 
An  important  input  for  this  calculation  is  the  magnetic  field  model. 
Ideally,  we  would  like  to  have  a  self-consistent,  time-dependent 
magnetic  field  model  (an  attempt  to  develop  such  a  model  based  on 
pressure  balance  between  field  and  particle  pressure  developed  numerical 
noise  and  was  postponed).  For  the  present  simulation  we  have  used  a 
superposition  of  the  Olson-Pf itzer  [1974]  analytic  model  and  a 
time-dependent  "substorm  current  loop"  that  simulates  the  effects  of  an 
induct -on  electric  field.  The  details  of  this  loop  will  be  discussed  in 
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Knowing  the  divergence  of  current  in  the  equatorial  plane,  we  solve 
equations  6  and  8  for  the  field-aligned  currents  into  the  ionosphere,  and 
equate  them  to  the  divergence  of  the  ionospheric  current  using  (7).  A 

semi-empirical  height-integrated  ionospheric  conductivity  model  has  been 

! 

developed,  details  of  which  are  presented  in  section  IV. C.  Equation  7  is 
an  elliptic  equation  in  two  dimensions  that  can  be  solved  numerically  for 
the  potential  V,  given  the  conductivity  tensor  and  j  y  ^  .  Our  boundary 
conditions  are  the  following: 

(i)  Zero  electric  current  across  the  equator  (this  condition  follows  from 
the  assumed  symmetry  between  the  hemispheres,  a  reasonable  assumption  for 
19  September,  which  is  near  equinox;  actually  the  condition  we  apply  is 
that  of  2ero  current  across  latitude  21°). 

(ii)  Specified  potential  V  on  the  polar  boundary.  Hereafter  (including 
papers  2  and  3)  we  use  the  term  "polar  boundary"  to  denote  the 
equatorward  boundary  of  the  region-1  Birkeland  current,  which  should  be 
distinguished  from  "polar  cap  boundary"  that  is  commonly  defined  by 
electric  field  reversals  (see  further  discussion  in  Section  IV .A). 

Because  of  the  irregularities  of  the  polar  boundary  we  actually  specified 
the  potential  on  a  circle  that  encompasses  the  polar  boundary.  The 
distribution  of  the  potential  around  this  curve  has  the  general  form 
suggested  by  Figure  1  of  Heelis  et  al.[1976];  the  magnitude  of  the 
potential  drop  was  estimated  from  real-time  observations  (see  section  IV. 
A). 

Given  the  ionospheric  potential  distribution  we  use  the 
magnetic-field  model  to  map  V  to  the  equatorial  plane.  In  the  present 


11 


simulation  we  neglect  the  component  of  E  that  is  parallel  to  B,  although 
we  believe  that  a  few  kV  potential  drop  over  limited  regions  would  not 
affect  our  results  substantially.  We  proceed  with  the  logic  loop  to 
calculate  magnetospheric  electric  fields.  The  total  electric  field  is 
the  sum  of  the  potential  electric  field  and  the  induction  field 
calculated  by  means  of  letting  the  equatorial  crossing  point  of  the  field 
line  vary  in  time.  This  motion  corresponds  to  E  x  B  drift  in  an 
induction  electric  field.  We  should  mention  that  this  is  only  one  of 
several  ways  to  introduce  the  induction  electric  field;  this  one  was 
chosen  for  its  simplicity  (see  equation  4  in  II. A). 

Given  the  potential  electric  field,  the  motion  of  equatorial 
crossing  points  due  to  induction,  and  the  magnetic-field  model,  the 
program  calculates  total  drift  velocities  (E  x  B  +  gradient  +  curvature) 
for  plasma-sheet  particles.  Specifically,  it  computes  the  motion  of  the 
inner  edge  of  each  species  k  of  the  plasma-sheet  electron-ion 
distribution,  namely  five  "energy  invariants"  for  electrons,  and  ten  for 
singly-charged  positive  ions. 

Loss  of  electrons  by  precipitation  is  included  in  the  present  model 
by  making  a  conventional  assumption,  namely  that  the  electrons  suffer 
strong  pitch-angle  scattering.  Under  these  conditions,  the  inner  edge  of 
the  electron  plasma  sheet  is  often  essentially  a  precipitation  boundary 
[Vasyliunas,  1968;  Kennel,  1969].  We  include  erosion  of  the  electron 
plasma  sheet's  inner  edge  in  an  approximate  way  that  involves  having  the 
computer  keep  track  of  two  boundaries  for  each  energy  invariant:  one 
boundary  where  25%  of  the  electrons  have  been  lost,  another  where  75X 
have  been  lost.  Proton  loss  has  been  ignored  in  this  first  set  of  runs. 

.. _ J 
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Given  the  velocities  of  different  components  of  the  inner  edge,  the 
program  advances  the  position  of  the  inner  edge  for  each  component  by  the 
amount  appropriate  to  one  time  step  fit  (30  seconds  in  these  runs)  by 
solving  equation  2.  This  closes  the  logical  loop  for  another  time  step, 
and  so  on.  i 

i 

In  the  actual  numerical  procedure,  the  program,  in  every'  time  step, 
reinterpolates  the  magnetic  field  model,  recalculates  Birkeland  currents 
for  the  new  particle  and  magnetic-field  configuration,  reads  the  observed 
electron  fluxes,  readjusts  Pedersen  and  Hall  conductivities,  readjusts 
the  polar  boundary  potential  drop,  re-solves  the  two-dimensional  elliptic 
equation  for  ionospheric  potentials  using  a  21  x  28  grid,  reinterpolates 
the  mapping  to  the  equatorial  plane,  calculates  corotation,  curvature  and 
gradient  drifts,  recomputes  boundary  velocities  and  moves  the  inner  edge 
of  various  components  of  the  plasma-sheet  (which  we  represent  by  400 
boundary  points). 

The  problem  is  further  complicated  by  the  fact  that  the  inner  edge 
of  the  plasma  sheet  is  often  rather  thin  (of  the  order  of  one  grid 
spacing).  Electric  fields  can  vary  by  large  factors  through  this  edge 
region  and  often  change  sign.  In  other  words,  electric  fields  generated 
by  one  part  of  the  inner  edge  strongly  affect  particle  motions  in  other 
parts.  To  accurately  model  this  sub-grid-scale  phenomenon,  we  have  had 
to  include  a  rather  intricate  self-correction  scheme,  which  substantially 
complicates  the  program. 


III.  DESCRIPTION  OF  THE  EVENT 


A.  Basic  description  of  the  event 

In  the  long  process  of  choosing  an  event  to  model,  we  had  to  scan 

many  data  bases^  for  the  years  1975-1976.  The  following  criteria  were 

* 

used  (some  of  them  somewhat  contradictory): 

a)  A  typical  and  clearly  defined  substorm  (with  one  major  expansion 
and  recovery)  so  that  our  simulation  model  would  be  as  general  as 
possible. 

b)  Long  enough  duration  so  that  we  could  have  several  passes  of 
polar-orbiting  satellites  during  the  course  of  the  substorm,  in  order  to 
maximize  available  observations  of  electric  fields,  magnetic  fields,  and 
auroral  particles.  (This  partly  contradicts  criterion  a). 

c)  Good  electric-field  data  from  several  passes  (before,  during  and 
after  the  substorm).  Several  passes  must  be  roughly  dawn-dusk  or 
dusk-dawn  to  allow  accurate  estimation  of  the  cross-polar-boundary 
potential. 

d)  As  many  good  DMSP  images  of  the  aurorae  as  possible. 

e)  Being  close  to  equinox  (to  minimize  difficulties  with  our 
untilted  magnetic  field  model). 

f)  North  America  being  at  nighttime  so  that  we  can  make  maximum  use 
of  North-American  ground  magnetometers.  This  implies  a  UT  restriction. 

g)  A  date  that  was  as  early  as  possible  so  that  the  data  could  be 
processed  and  supplied  to  us. 

h)  Availability  of  other  magnetospheric  and  ionospheric  data  (to  be 
used  in  future  detailed  analysis  of  the  event). 


i)  Positioning  of  a  satellite  near  the  inner  edge  of  the  plasma 
sheet. 

The  substorm  we  finally  chose  satisfied  criteria  b,  c,  e,  f,  g  and 
partly  satisfied  criteria  a,  d  and  h. 

Figure  3  (top  panel)  shows  the  Fort  Churchill  magnetogram  for  19 
September  1976  as  a  function  of  Universal  Time  (Greenwich  Mean  Time). 

The  event  we  have  chosen  to  model  is  the  one  with  an  onset  at  1000  UT  and 
a  duration  of  about  5  hours;  at  onset  Fort  Churchill  was  at  1455  MLT.  A 
small  disturbance  was  observed  a  few  hours  earlier  (around  0600  UT); 
unfortunately,  two  of  the  S3-2  satellite  passes  occurred  during  that  time 
span. 

Further  investigation  of  auroral  and  mid-latitude  magnetograms 
reveals  that  there  was  more  than  one  substorm  in  the  time  span  of  our 
simulation.  Virtually  all  nightside  stations  show  a  large  substorm  with 
onset  at  about  1000  UT  or  shortly  thereafter,  but  some  stations  (e.g. 
College  and  Pacific  Ocean  stations)  show  a  smaller  localized  expansion  at 
1150  UT.  We  cannot  even  rule  out  the  possibility  of  additional  minor  and 
localized  disturbance.  Our  method  of  computing  conductivities  and  the 
boundary  potential  drop,  based  on  discrete  satellite  passes,  cannot 
resolve  these  small  recoveries  and  expansions.  Therefore  we  computed  the 
response  of  the  inner  magnetosphere  to  one  substorm  with  a  long,  slow 
recovery  although  it  appears  that  the  actual  situation  was  considerably 
more  complex. 

B.  Satellite  data 

The  primary  satellite  used  for  the  ionospheric  electric  field  was 
the  Air  Force's  S3-2  satellite.  This  satellite  is  a  spin-stabilized 
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polar  orbiter.  Inclination  97°.  It  has  two  dipole  electric-field  probes, 
one  in  the  spin  plane,  which  is  nominally  the  orbital  plane,  and  one  a- 
long  the  spin  axis.  It  also  has  a  tri-axlal,  fluxgate  magnetometer  and 
an  energetic  electron  detector.  Data  were  stored  on  an  onboard  tape  re¬ 
corder.  Coverage  of  data  was  limited  by  the  relatively  few  tracking  sta¬ 
tions  available  to  transmit  the  data. 

Figures  4-6  present  electric  fields,  magnetometer  and  electron-ener¬ 
gy-flux  data  from  the  S3-2  satellite  for  3  passes  during  the  event.  (For 
further  discussion  of  the  processing  of  the  data,  see  Burke  et  al  [1979]). 

The  component  of  electric  field  in  the  forward  direction  of  the  sa¬ 
tellite  and  approximately  perpendicular  to  the  magnetic  field  is  shown  in 
the  bottom  panels  of  Figures  4-6,  The  conventional  polar  cap  potential 
drop  (boxes  in  Fig. 3,  bottom)  was  obtained  by  integrating  E  •  dl  across  the 
polar-cap  (from  reversal  to  reversal).  As  is  obvious  from  Figures  4-6, 
the  polar  cap  boundary  is  sometimes  hard  to  determine;  the  most  equator- 
ward  reversal  was  generally  used  (arrows  in  Figures  4-6).  The  middle  pan¬ 
el  of  Figures  4-6  show  the  magnetic  variations,  AB  along  the  spin  axis  of 
the  satellite,  which  was  almost  exactly  east-west.  It  is  apparent  that 
region-l-type  currents  often  extended  well  into  the  sunward  convection  re¬ 
gion  (denoted  by  vertical  lines,  e.g.  1134  UT  in  Figure  6).  Thus  the  "po¬ 
lar  boundary  potential  drop",  from  the  viewpoint  of  our  model  (lines  with 
error  bars  in  Figure  3b),  is  often  significantly  smaller  than  the  conven¬ 
tional  polar  cap  potential  drop. 

Figures  4-6  are  significant  in  that  they  are  among  the  first  satellite 
data  across  the  polar  cap  showing  simultaneous  aB  perturbations  and  electric 
field  data.  These  data  reveal  some  Interesting  features:  there  is  a  general 
correlation  between  the  magnitude  of  the  transverse  magnetic  field  perturba¬ 
tion  aB  and  the  forward  comDonent  of  elertrir 
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field.  Simultaneous  measurements  of  particles  and  magnetic  fields  on  the 
ISIS  2  satellite  during  quiet  times  were  found  to  be  correlated  [Burrows 
et  al.,  1976;  Klumpar  et  al.,  1976;  and  McDiarmid  et  al.t  1977,  1978]. 
McDiarmid  et  al.  [1978]  further  compared  magnetic  perturbation  profiles 
with  average  electric  field  profiles  reported  by  Heppner  [1972]  and 
Gurnett  [1972]  and  found  the  shapes  to  be  very  similar,  suggesting  that 
the  magnetic  field  tilts  in  the  direction  of  the  convection.  However,  in 
two  of  the  passes  presented  here  (Figures  5  and  6)  regions  of  strong 
sunward  convection  are  observed  poleward  of  the  strongest  region-1 
currents,  suggesting  that  despite  the  general  similarity,  one  cannot 
predict,  at  least  during  substorms,  the  plasma  convection  based  entirely 
on  AB,  as  suggested  for  steady-state  conditions  [McDiarmid  et  al.,  1978]. 

The  auroral  region  electric  field  and  magnetic  perturbation  data  in 
Figures  4-6  were  used  for  comparison  with  theoretical  results.  Detailed 
discussion  and  comparison  with  theoretical  predictions  are  presented  in 
the  following  paper. 

Also  shown  on  Figure  5  are  features  of  the  Lockheed  Ion  Mass 
Spectrometer  data  from  a  nearly  simultaneous  S3-3  pass  (R.  D.  Sharp, 
personal  communication,  1977).  The  S3 -3  pass  was  also  dawn-dusk,  from  an 
invariant  latitude  of  61.6°,  magnetic  local  time  of  1805  at  1005  UT,  to 
an  invariant  latitude  of  54°,  magnetic  local  time  of  0730  at  UOOUT.  The 
S3-3  polar  cap  pass  covered  a  much  longer  time  span  because  of  its  higher 
altitude  (which  dropped  from  7806  km  to  2704  km  during  this  time 
period).  Despite  minor  differences  in  MLT  and  UT  between  S3-2  and  S3-3, 
the  electron  fluxes  observed  at  corresponding  invariant  latitudes  were 
remarkably  similar.  In  particular,  precipitating  electron  bursts 


17 


resembling  "Inverted  Vs"  were  observed  near  68°,  69°,  74®,  and  77°  In¬ 
variant  from  both  satellites.  In  the  S3-2  data,  these  bursts  were  associated 
with  upward  field-aligned  currents;  in  the  S3-3  data-*  these  bursts  were 
associated  with  upward  flowing  ions.  Thus  it  seems  reasonable  to  infer  a 
connection  between  upward  flowing  iops  and  upward  Birkeland  currents.  It 

i 

is  reassuring  to  note  that  no  upward'  flowing  ions  were  observed  in  our 
modeling  regions  (not  even  in  the  upward  current  region  on  the  dawns ide), 
since  upmrd  flowing  ions  may  be  associated  with  parallel  electric  fields 
[Shelley  et  al..  1976;  Mizera  and  Fennell.  1977;  Ghielmetti  et  al., 

1978],  and  we  have  assumed  no  parallel  electric  fields  in  our  modeling 
region. 


I 
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IV.  INPUT 

A.  Polar  boundary  potential  distribution 

The  polar-boundary  potential  drop  is,  in  our  view,  the  most 
important  single  input  parameter  to  our  model.  The  potential  drop  is  our 
monitor  of  the  solar  wind  and  boundary  layer's  influence  on  the  inner 
magnetosphere.  It  should  be  emphasized  that  we  have  taken  our  polar 
boundary  at  the  equatorward  edge  of  the  region-1  current,  in  contrast 
with  the  more  conventional  polar  cap  defined  by  the  electric  field 
reversal  (see  Figure  3,  bottom). 

Measurements  shown  in  Figure  3  (bottom)  indicate  that  the  potential 
drop  is  relatively  constant  prior  to  onset,  and  increases  during  the 
growth  and  main  phases  of  the  substorm.  (Such  an  increase  was  also 
suggested  by  Mozer  [1973],  based  on  balloon  measurements  ).  However,  as 
a  result  of  the  data  points  being  so  far  apart  it  is  difficult  to  derive 
a  "best  fit"  curve.  The  curve  we  have  chosen  is  shown  in  solid  line  in 
Figure  3  (bottom  panel).  The  potential  drop  was  kept  constant  at  20  kV 
utill  0900  GMT.  Then  we  increased  it  linearly  to  reach  a  peak  of  80  kV 
at  1050,  after  which  it  stayed  constant  throughout  the  long  recovery 
phase  of  our  model  substorm. 

As  a  result  of  the  uncertainty  in  the  time  dependence  of  the  polar 
boundary  potential  drop,  we  ran  our  program  again  with  a  different 
assumption  about  the  potential  drop  (dotted  line),  with  a  peak  value  of 
140  kV.  Detailed  comparison  with  data  shows  the  smaller  potential  drop 
to  be  more  realistic. 
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In  addition,  the  following  difficulties  were  encountered  trying  to 
derive  a  potential  drop  as  a  function  of  time  during  the  event: 

(1)  Although  most  of  the  passes  were  mainly  dawn-dusk,  some  only 
"skimmed"  the  polar  cap  boundary  {in  particular  the  first  two)  giving 
rise  to  an  error  (probably  an  underestimate)  in  the  resulting  potential 
drop. 

(ii)  We  have  no  potential -drop  measurements  between  0640  and  0940  GMT. 
Combined  with  the  fact  that  the  0450  and  0620  GMT  passes  were  during  the 
small  disturbance  and  were  therefore  ignored,  we  had  virtually  no  data 
between  0450  and  1000  (onset).  As  a  result  we  believe  that  we  may  have 
underestimated  the  length  of  the  growth  phase  for  this  run,  represented 
by  the  solid  curve  in  Figure  3  from  0900  to  1000  GMT. 

Given  the  potential  drop,  we  apply  a  simple  analytic  function  to 
estimate  the  potential  distribution  on  the  polar  boundary.  The 
electric  equi potential  curves  were  compressed  across  the  dayside  in 
accordance  with  observations  [Heelis  et  a!.,  1976].  The  polar 
boundary  was  assumed  to  be  an  offset  circle  [e.g. ,  Feldstein,  1973;  Meng 
et  al.,  1977].  We  have  chosen  an  offset  of  two  degrees  toward  midnight, 
away  from  the  geomagnetic  pole.  The  results  of  the  first  simulation 
suggest  that  our  results  can  be  improved  considerably  by  a  more 
judicious  choice  for  the  fit  to  the  polar  boundary  potential  distribution. 
B.  Magnetic  field  model 

For  these  simulation  runs  we  used  the  analytic  version  of  the 
Olson-Pf itzer  [1974]  magnetic  field  model.  The  model  appears  basically 
adequate  for  our  purposes,  although  It  has  the  disadvantage  that 
typically  |V*||  -  Sy/R^  on  ring-current  field  lines.  To  the  Olson- 
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Pfitzer  analytic  model,  we  added  the  magnetic  field  of  a  substorm  current 
loop  (see  also  section  II),  which  is  used  to  estimate  the  time-dependent 
magnetic  field  (and  the  resulting  induction  electric  field)  (Fig.  7). 

The  loop  was  first  suggested  by  Atkinson  [1967]  and  used  by  McPherron  et 
aK ! C 1973]  and  many  others.  Our  version  of  the  loop  has  an  equatorial 
eastward  current  across  the  tail  at  50  (representing  interruption  of 
the  tail  current  during  a  substorm),  Birkeland  currents  down  to  the 
northern  and  southern  hemisphere  ionospheres,  westward  electrojet,  and 
another  Birkeland  current  from  the  ionosphere  out  to  the  equatorial 
plane.  For  this  first  run,  we  chose  the  maximum  current  in  the  loop  such 
that  it  would  approximately  reproduce  the  ground  magnetic  variation 
(positive  bay)  observed  at  Boulder.  The  loop  was  turned  on  at  onset  and 
reached  its  maximun  value  at  1040  UT  (peak  of  the  substorm)  after  which 
it  gradually  declined  through  the  long  model  recovery  phase. 

Figure  8  shows  contours  of  constant  magnetic  field  in  the  equatorial 
plane  for  0900  UT  (  an  hour  before  onset).  Although  symmetric  and 
dipole-like  near  earth,  the  contours  deviate  from  this  symmetry  for 
larger  L  values.  The  contours  for  L  <  6  will  be  essentially  the  same 
even  later  in  the  event,  when  the  substorm  loop  is  at  full  strength.  For 
larger  L  values,  the  contours  near  midnight  are  pulled  antisunward, 
representing  return  to  a  more  dipole-like  structure  near  midnight. 

C.  Conductivity  Model 

The  conductivity  model  is  an  important  input  to  the  equation  for 
conservation  of  ionospheric  current  (equation  4  in  section  II),  which  we 
solve  numerically  to  find  the  potential  distribution  in  the  ionosphere. 
The  height-integrated  Pedersen  and  Hall  conductivities  include  base-level 
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and  time-dependent  terms.  The  base- level  term  includes  day-night 
asymmetry  and  solar-zenith-angle  dependence.  Nightside  midlatitude 
values  are  estimated  from  Rowe  and  Mathews,  [1973]  and  roughly  agree  with 
Harper  and  Walker  [1977];  dayside  conductivities  are  consistent  with 
Harper  [1977]. 

To  estimate  the  spatial  dependence  of  the  auroral  enhancement  in  an 
accurate  way,  for  the  whole  auroral  zone,  one  would  need  global 
measurements  of  electron  density  or  global  measurements  of  the  flux  of 
auroral  electrons  incident  on  the  ionosphere.  Such  measurements  are,  of 
course,  not  available,  and  one  must  either  approximate  or  extrapolate 
available  data.  Our  present  conductivity  model  is  a  first  attempt  at 
including  these  enhancements  in  a  time-dependent  way. 

The  primary  method  we  have  used  to  estimate  height-integrated 
conductivities  for  the  19  September  event  involves  an  empirical  formula 
that  we  derived,  relating  height-integrated  Pedersen  and  Hall 
conductivities  to  auroral  electron  energy  flux  and  mean  electron  energy. 
M.H.  Rees  and  his  collaborators  [Jones  and  Rees,  1973,  and  Rees,  private 
communication]  have  calculated  electron-density  profiles  in  the 
night-time  auroral  ionosphere  for  various  energy  fluxes  and  mean 
energies,  and  for  assumed  neutral -atmosphere  models.  Using  standard 
formulas  for  Pedersen  and  Hall  conductivities  [Rishbeth  and  Garriott, 
1969],  we  deduced  height  profiles  of  Hall  and  Pedersen  conductivities  for 
three  of  Rees'  models.  We  then  numerically  integrated  these  over  height 
and  fitted  the  results  to  simple  power  laws,  with  the  following  results: 
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4Eh  =  0.55  [Average  Electron  Energy  f •«  __  ^  ,10) 

Energy  flux  and  average  electron  energy  were  computed  from  S3-2 
electron  data,  using  32  electron  channels  that  cover  the  range  .08-17  keV. 

For  our  substorm  simulation  we  used  flux  measurements  from  four 
passes  on  19  September  1976:  0400,  1000,  1150,  and  1500  GMT.  For  each 

pass  we  determined  the  flux  as  a  function  of  invariant  latitude. 

We  have  constructed  conductivity  models  with  two  different  degrees 
of  latitudinal  smoothing  at  the  inner  edge  of  the  plasma  sheet  (Table 
2).  The  smoothing  is  necessary  because  of  the  inability  of  the  current 
conservation  difference  equation  to  cope  with  very  sharp  conductivity 
jumps.  For  a  detailed  numerical  description  of  the  conductivity  model, 
see  section  V. 

Figure  9  demonstrates  the  latitudinal  dependence  of  our  conductivity 
model  for  these  "less  smooth"  conductivities.  We  plot,  for  a  given  time 
(1150  UT),  the  Pedersen  and  Hall  conductivities  as  a  function  of  local 
time  for  4  different  latitudes.  Note  the  conductivity  peaks  at  local 
noon  except  for  the  70°  profile  when  the  auroral  enhancement  causes  the 
conductivity  to  peak  on  the  night  side. 

Figure  10  shows  the  universal  time  dependence  of  the  conductivity 
model.  We  plot  the  Pedersen  conductivity  profile  at  auroral  zone 
latitude  (70°  )  for  various  phases  of  the  substorm.  The  conductivity 
peak  coincides  with  the  peak  in  the  cross-polar-cap  potential  drop  (at 
1150  UT). 
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D.  Initial  hot-particle  distribution  in  the  magnetosphere 

The  initial  plasma  configuration  that  we  assumed  is  somewhat 
arbitrary,  since  no  magnetospheric  hot-plasma  data  were  available.  The 
inner  edge  of  the  plasma  sheet  was  assumed  to  have  a  finite  thickness 
(about  0.7  )  and  was  located  at  L  «10.  An  isotropic  pitch  angle 

_3 

distribution  was  assumed  with  nfi  =  0.9  cm  ,  kTg  =  1.5  keV  and  kT^  = 

4.5  keV  at  L  =10.  Unlike  the  other  input  parameters,  these  values  were 

inserted  once  as  initial  conditions  and  were  calculated  at  later  times 

2/3 

with  the  assumption  that  E^  (/  ds/B)  '  serves  as  an  adiabatic 
invariant  (  =  kinetic  energy  of  particles  of  type  k,  /ds/B  =  flux 

tube  vo 1  une). 


n 
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V.  NUMERICAL  METHOD 

The  Earth's  magnetosphere  is,  of  course,  an  intrinsically  three- 
dimensional  system  that  cannot  be  realistically  approximated  as  a 
two-dimensional  system,  uniform  in  the  third  dimension.  One  approach 
to  the  problem  is  therefore  to  attempt  to  solve  the  simultaneous  time- 
dependent  three-dimensional  MHD  equations.  One  group  is  attempting  to  do 
this  [Lebouf  et  al . ,  1980].  Our  approach  has  been  to 
reduce  the  problem  to  a  treatable  two-dimensional  computation  by  using 
the  fact  that  the  magnetic-field  configuration  in  the  inner  magnetosphere 
can  be  estimated  relatively  accurately,  and  the  field  lines  are  relatively 
good  conductors.  In  fact,  we  solve  two  2-dimensional  sets  of  equations 
(in  the  ionosphere  and  in  the  Earth's  equatorial  plane);  the  two  regions  are 
coupled  via  the  geomagnetic  field. 

The  system  of  equations  (1  -  8)  in  section  II. A  describes  the  motion 
of  the  plasma  sheet  during  disturbances  and  the  resulting  ionospheric 
currents  and  electric  fields.  Solution  of  the  system  of  equations  involves 
an  intricate  and  elaborate  numerical  analysis.  We  therefore  will  describe 
only  some  of  the  most  important  aspects  of  the  numerical  treatments. 

A.  Coordinate  System 

In  the  ionosphere  we  have  chosen  an  orthogonal  2-dimensional  coordinate 
system  (c.ip)  similar  to  the  one  described  in  Jaggi  and  Wolf 
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[1973].  The  center  of  the  coordinate  system,  however.  Is  displaced  2 
degrees  toward  midnight  from  the  geomagnetic  north  pole  (so  that  our 
polar  boundary,  which  is  shown  in  Figure  lb  and  is  chosen  to  be 
independent  of  time,  follows  a  contour  of  £  *  constant).  A  length 
increment  ds  in  the  ionosphere  can  thus  be  described  as 

ds2  =  a2  d£2  +  B2d^2  (11) 

where  a(£)  and  8(£)  are  some  scalar  functions.  The  resulting  grid  system 
(i,j)  consists  of  18  x  28  grid  points  covering  all  magnetic  local  times. 
The  points  are  evenly  spread  in  j  with  j  =  3  or  31  at  local  noon,  j  =  10 
at  dusk  and  j  =  17  at  midnight.  The  grid  covers  approximately  the  range 
20°<  0  <  72°  in  colatitude.  The  latitudinal  spacing  of  grid  points 
a(£)A£  varies;  in  the  main  convection  zone,  a(£)A£  is  equivalent  to  about 
1.7°,  while  at  mid-  and  low  latitudes  it  varies  from  *  1.8°  to  =  8.7°. 

The  local-time  grid  spacing  8A\|»  is  equivalent  to  about  0.86  hours. 

The  equatorial  grid  system  is  obtained  by  tracing  magnetic  field 
lines  from  ionospheric  grid  points  out  to  the  equatorial  plane.  This 
results  in  a  nonorthogonal  grid  for  the  equatorial  plane.  The  equatorial 
extent  of  our  grid  system  is  roughly  1.03  <  L  <  10  and  is  somewhat 
time-dependent  (since  our  magnetic  field  model  is  also  time-dependent). 

B.  Inputs 

Aside  from  the  grid  system  various  input  parameters  are  needed  to 
solve  the  time  dependent  system  of  difference  equations  that  represent 
equations  1-8  in  section  II. A. 

(i)  Time  Independent  Inputs 

For  an  isotropic  pitch  angle  distribution,  the  distribution  function 


Is  given  by 


2.6 


f(v)  « 


-  mv  /2kT 


(2tt*cT  /m) 


3/2 
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where  n  is  the  nunber  density. 

Using  X  (the  "energy  invariant")  and  n  (the  "density  invariant"),  as 
defined  in  section  II. A.,  one  can  derive 


2nT  A*  e"A/Vr 

dn/dx  =  — L - -372_  (13) 

*  X<T 

where  is  the  total  nunber  of  particles  per  unit  magnetic  flux  in  the 
flux  tube  and  X^  =  <T(/ds/B)273.  We  keep  track  of  15  distinct  "species". 
For  the  electrons. 


=  2579.3  sinh  [0.1925  (k-6)]  1  <  k  <  5  (14) 

\  ~  Xk-5  6  <  k  <  10 

and  for  the  ions 

Vqk  =  317'°  Sinh  t0'4  11  <  k  <  21 

where  X  is  in  units  of  eV  Rg273y"273J  k  is  a  running  index  and  is  the 
particle's  charge.  We  use  two  k  values  for  each  electron  X  value,  as 
part  of  our  procedure  for  including  electron  loss  by  precipitation.  (See 
subsection  C,iii  below.)  The  k  =  11  species  represents  cold  plasma.  The 
resulting  densities  are,  for  electrons, 

nk  =  %  (-^r— )  (-,—)*  cosh[0.1925(k-6  )]  exp  [-  y^-]  l<k<5 

K  1  ATe  ATe  ATe 


nk  =  nk-5 


6  <  k  <  10 


(15) 
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and  for  ions 


£43.Jj  (  Xk  )*  cosh  [0.4  (k-11)]  exp  [-“7-]  11<  k<2 


nk  "  nT  rXTi 


The  constant  parameters  are: 

21  -1 
»j-  =1 .7x  10  weber 


XTe=  721  eVy“2/^E2/3;  X^.  =  2163  eVy‘2/3RE2/3 
(ii)  Time  Dependent  Inputs 
a)  Magnetic  Field 

The  magnetic-field  model  at  a  given  time  t  is  described  numerically 
by  three  matrices:  ;  R1J;  P1J.  For  a  given  i,  j  ionospheric  grid 

point,  these  matrices  give  the  Bz  and  the  polar  coordinates, 
respectively,  of  the  corresponding  point  in  the  equatorial  plane  (the 
mapping  itself  is  done  via  the  magnetic  field).  For  the  period  before 
substorm  onset,  the  matrices  were  calculated  using  the  analytic  version 
of  the  Olson-Pfitzer  model  [1974]  to  trace  each  grid  point  from  its 
ionospheric  origin  to  its  equatorial  crossing.  For  the  peak  of  the 
substorm  this  field-line  tracing  procedure  was  repeated  using  a  modified 
magnetic  field  that  included  the  effect  of  the  substorm  current  loop  (see 
section  IV. B).  For  other  times,  we  interpolated  between  the 
"quiet-time"  and  the  "substorm-peak"  values  of  the  three  matrices.  The 
flux  tube  volume  was  approximated  using  the  formula 


(/#)lj  =  [Baij/Rij  +  20  (R^')"1'3]"1 


b)  Conductivity 
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approximation  of  layer  conductivity  where 

E«  -  Vs1"21 

=  "  *>£  = 
v  =  r 

n  p 

where  E^,  E^,  E  ^  and  E^  are, the  components  of  the  conductivity  tensor 
in  equation  7,  Ep  and  E^  are  height-integrated  Pedersen  and  Hall 
conductivities,  and  I  is  the  magnetic  dip  angle.  (This  approximation 
holds  if  we  stay  away  from  the  geomagnetic  equator).  To  calculate  the 
height-integrated  conductivities  in  mhos  at  grid  point  (i,j)  at  time  t, 
we  use  the  following  formulas 


(0  - 1 

K  hemispheres 

(O  =  l 

hemi spheres 


|  12.5  S[^cosxij,-u(cosxij)]+0.45  +  AEplj(t)  J  (18) 
|  18.75  S[^cosx^*  -u(cosxlj  )]-*&#*  +  AEH1J(t)|  (19) 


In  each  formula,  the  first  term  represents  the  time- independent  sunlight 
effect,  X1J  =  solar-zenith  angle  at  (i,j),  computed  assuming  that  the 
earth's  magnetic  dipole  is  perpendicular  to  the  earth-sun  line;  u(x)  =  1 
or  0  depending  on  whether  x  >  0  or  x  <  0;  and  S  represents  a  numerical 
function  that  tries  to  smooth  slightly  the  sharp  jumps  at  dawn  and  dusk. 
The  second  terms  in  (18)  and  (19)  give  base-level  nightside  values.  The 
third  terms  represent  time-dependent  auroral  enhancements;  they  were 
computed  in  a  complicated  way  from  (9)  and  (10),  which  give  the 
enhancements  below  the  S3-2  satellite,  by  linearly  interpolating  in 
universal  time  to  get  values  for  times  t  between  satellite  passes, 
extrapolating  smoothly  in  local  time  assuming  twice  as  great  enhancement 
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at  midnight  as  at  noon,  converting  from  actual  GM  latitude  to  an 
effective  model  latitude,  and  smoothing  latitudinal  fine  structure. 

c)  Polar  boundary  potential  distribution 

S3 -2  satellite  measurements  were  used  to  estimate  the  dawn-dusk 
potential  drop  across  the  polar  cap  (see,  e.g..  Figure  4-6).  The  time 
dependence  was  approximated,  again,  by  linear  interpolation  (Figure  3b). 
The  longitudinal  dependence  was  chosen  so  as  to  concentrate  the  potential 
drop  across  a  smaller  distance  on  the  dayside  (the  "throat"  region 
[Heel is  et  al.,  1976])  than  on  the  nightside. 

The  boundary  potential  is  given  by 

Vb  =  Vosin  near  noon  1  */4  or  >  7tt/4 ) 

=  VQ  for  ir/4  <  ifi  <  ir/2 

Vb  =  VQsim|>  for  the  night  sector  (ir/2  <.  <  3n/2)  (20) 

Vb  =  -VQ  for  3tt/2  <  <|>  <  7tt/4 

where  Vq  *  1/2  x  (cross  polar  boundary  potential  drop). 

C.  Calculations  of  electric  fields  and  plasma  motion 

Our  general  approach  to  the  solution  of  our  problem  is  schematically 
shown  in  Figure  2  and  discussed  in  section  II. B.  In  order  to  complete  a 
time  step  in  our  logic  loop  we  have  to  solve,  simultaneously,  equations 
1-8.  We  split  this  task  into  3  parts;  calculation  of  the  potential 
distribution  V  in  the  ionosphere  (equation  7);  calculation  of  the 
electric  field  at  an  arbitrary  point  and  calculation  of  the  plasma  sheet 
motion  (equation  2). 

(i).  Calculation  of  the  potential  matrix  V 

The  electrostatic  potential  V  is  computed  using  the  conservation  law 
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for  ionospheric  currents  (equation  7).  Once  we  compute  the  horizontal 
conductivity  and  the  Birkeland  currents  we  approximate  equation  7 
numerically  by 


V 


(21) 


where  (i£.j£)  =  0+1. j).  (i»j+l)  and  (i,j-l),  respectively,  for 

1=1,  2,  3,  4.  The  factors  Cj  -  are  the  same  as  given  in  Jaggi  and 
Wolf  [1973],  and  C^J  is  given  by 


ij  .  Vf  1'j 
i  "  iL5k 
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(22) 


The  present  case  is  different  from  the  case  presented  in  Jaggi  and  Wolf 
[1973]  in  that  many  “species"  are  now  considered  (running  index  k  over  21 
"species").  The  term  AV^  =  amount  by  which  the  potential  at  (i£,j£)  is 
changed  due  to  presence  of  discontinuity  k,  and  (i0»jQ)  =  (i,j).  As 
discussed  later,  the  sun  includes  only  layers  that  cut  the  "cross", 
namely  the  lines  connecting  (i,j)  to  the  (i£*j£)»  *  s  1,..,4.  We  find 
that 


AV 


-  I  I 

tk  "  2  L*k 


_ _  ] 

cos2;  +  fy,s1nzc 


(23) 


plus  a  constant  that  does  not  depend  on  *  and  eventually  cancels  out. 
Here  =  perpendicular  distance  from  inner  edge  k  to  grid  point 
(i^.j^),  j  ij  ^  =  Birkeland  current  coming  down  into  the  ionosphere  along 
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the  inner  edge  for  species  k,  per  unit  length  along  the  inner  edge,  and 
C*  angle  between  the  inner  edge  and  the  i  =  constant  grid  line.  (Equation 
23  is  based  on  the  assumption  that  the  conductivity,  and  the  electric 
field  due  to- a  layer, is  approximately  constant  over  the  order  of  a  grid  spacing). 
Substituting  equation  23  into  equation  22  gives 
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One  problem  here  is  that  in  computing  the  AV's,  we  neglected  conductivity 
inhomogeneities.  It  turns  out  to  be  better  in  the  formula  for  C5  to  also 
neglect  the  conductivity-gradient  terms  in  the  C^s,  namely  to  make  the 
following  replacements,  using  the  notation  of  the  Jaggi-Wolf  paper: 


C  =  —  V 
5  4  l  r 

k  i 


ilLL 


rr  cosZc+I 
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ij 


(25) 


A^r 


The  quantity  in  curly  brackets  can  be  shown  to  be  exactly  zero  for  any 
edge  that  fails  to  cut  the  cross,  assuming  a  straight  boundary  and  equal  grid 
spacing,  and  thus  we  can  self-consistently  consider  the  sum  to  be  restricted 
to  species  values  whose  inner  edges  cut  the  cross. 

The  sum  in  (25)  is  a  numerical  approximation  to  an  integral  over 
energy  invariant  X.  Ideally,  we  would  have  enough  species  k  that  many 
k-values  would  break  the  cross  and  contribute  to  the  sum  for  any  (i,j). 

If  the  spacing  between  inner  edges  of  different  k  happens  to  be  large, 
then  the  sum  for  a  given  (i,j)  may  be  a  poor  approximation  to  the 
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corresponding  integral  over  A  for  that  (i,j).  To  take  account  of  this 
problem,  we  define  I  layer  subsPec'es  between  eachpairof  successive 
k-values  (species)  and  attribute  to  each  subspecies?  inner  edge  an 
interpolated  position  and  a  reduced  and  interpolated  subspecies  jjj. 
Subspecies  are  labeled  by  the  index  k '  •  The  sun  in  (25)  then  becomes  a 
sum  over  subspecies  k*.  The  choice  of  I-|ayer  is  discussed  in  subsection  E. 

The  difference  equation  (21 ),  representing  equation  7  in  section  II. A, 
is  solved  by  a  straightforward  iterative  technique  (see  also  Oaggi  and 
Wolf  [1973]).  For  a  given  time  t,  we  derive  the  boundary  potential  VJ  on  the 
polar  boundary,  as  described  in  subsection  B.ii.c.  Given  the  boundary 
condition  =  0  at  the  equatorward  edge  (lat.  «  21°),  given  a  set  of 
values  V1J  as  a  "guess"  and  computing  the  coefficients  Cj1J  -  we 

obtain  an  improved  set  of  values  (V')1J  as  follows: 


6  Vij  =W[C1ij  Vi+1’jK2ij  V1"1^  +  C3ij  VUj+l 
+  C4ij  -  Vij] 

( V  ' ) 1  J  =  V1J  +  6Vij 


(26) 

(27) 


where  W  is  a  weight  factor  (0  <  M  s  1).  This  procedure  is  repeated  many 
times  until 

l  |6Vlj|2<e  (28) 

Provided  the  initial  guess  is  good,  criterion  (28)  is  satisfied  within 
'50  iterative  steps.  We  should  comment  that  this  simple  iteration 
scheme  is  efficient  because  we  have  a  very  good  initial  guess  for  V  each 
time  step,  based  on  V  for  the  last  two  time  steps. 

(ii)  Calculation  of  the  electric  field 

The  simplest  and  most  straightforward  manner  of  computing  the 
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electric  field  at  an  arbitrary  point  is,  of  course,  to  use  some  form  of 
interpolation  based  on  the  potentials  V*J  at  the  grid  points.  This 
method,  however,  will  not  suffice  in  the  region  of  the  "inner  edges"  due 
to  the  large  electric-field  discontinuities  at  some  inner  edges.  We  thus 
developed  a  more  intricate  and  complex  method  of  computing  the  electric 
field. 

To  calculate  motions  of  the  inner  edges,  we  must  do  accurate  computations 
of  electric  fields  at  noninteger  i-values,  which  we  label  I.  We  start  from 
an  interpolation  formula,  then  correct  it  for  sharp  inner  edges.  To  compute 
the  two  components  of  W  at  the  point  (I,j),  where  i  <  I  <  i+i,  we  use  the 
following  formulae: 


fj(l.j)  =  \  [Vi+1>j  -  +  (I-i)[Vi+1,i  -  2V1 +  V1-'’5] 

+  lt  \  [AVk,i+1’j  -  AVk.1_1*j]  -  (I-i)[AVk,1+1*J‘  -  2AV k , 1  ’ ^  ♦  AVk,i-1’j 

k  ...  (29) 

=  |cv1,j+1  -  V1”*'1]  +  (I-i)  [V1+1,j+1  +  V1,J  -  V1,J+1  -  V1+1,j] 

+  -  ?  UVk!’J+1  -  AVk.i»j’1]  -  (I-i)[AVk,i+1’j+1  +  AVk,i,j 

-AVk,i,j+1  -  AVk.i+1,j]}  (30) 


Here  and  later  in  the  paper,  we  use  I  and  J  for  non-integer  values  of  the 
grid  indices  i  and  j.  The  firstlinesof  (29)  and  (30)  give  a  straight¬ 
forward  interpolation  formula.  In  each  second  line,  the  first  terms 
give  the  appropriate  component  of  the  potential  gradient  generated  by 
each  inner  edge  k'.  The  rest  of  each  second  line  and  third  line  subtracts  off  the 
interpolation-formula  result  for  the  potential  gradient.  Thus  each 
curly-bracketed  term  gives  the  difference  between  the  potential  gradient 

at  (I,j)  due  to  inner  edge  k1  and  the  approximation  to  that  potential 

i .  i 
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by  (23),  is  proportional  to  the  distance  from  grid  point  (i,j)  to  inner 
edge  k'  and  to  the  Birkeland  current  per  unit  length  along  inner  edge  k* 
as  it  crosses  grid  line  j.  The  curly-bracketed  correction  terms  in  (29) 
and  (30)  are  zero  if  the  point  (I.j)  and  all  the  relevant  grid  points  lie 
on  the  same  side  of  inner  edge  k'.  , 

i 

(iii).  Plasma  Motion 

The  inner  edge  of  a  given  species  k  is  specified  by  a^(t),  where  j 
is  our  local  time  index,  which  varies  between  j  =  3  and  j  =  30.  That  is, 
we  use  the  point  where  the  inner  edge  k  crosses  the  j  =  integer  line 
(these  are  lines  of  =  constant  longitude).  In  most  cases,  the  inner 
edges  tend  to  run  almost  parallel  to  lines  of  constant  latitude.  Our 
task  is  to  compute  3a^(t)/3t. 

The  basic  formula  for  drift  velocity  of  particles  E  x  B  drifting  in 
a  potential  electric  field  and  also  gradient  and  curvature  drifting 
(assuming  isotropic  pitch  angles)  is  the  following: 

vD  =  B  x  Weff/B2  (31) 

where  Vg^^includes  the  convection  potential,  the  corotation  potential  and 
the  effects  of  gradient  and  curvature  drifts  (see  equations  2  and  3  in 
section  II. A).  If  we  now  apply  this  formula  to  the  ionosphere  and  use  it 
to  calculate  the  rates  of  change  of  i  and  j  coordinates  due  to  drift,  we 
find 
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where  Bjr  =  radial  component  of  B  at  the  ionosphere.  In  deriving  these 
formulas,  we  have  had  to  take  into  account  the  fact  that  the  magnetic 
field  is  not  exactly  perpendicular  to  our  ionospheric  grid,  which  lies  on 
a  spherical  shell.  The  rate  of  change  of  the  i-value  of  the  intersection 
of  layer  k  with  grid  line  g  is  given  by 

9ak  _  dl_  _  dj  aak 
3t  dt  dt  9J 

3Vof-  3Voff  3a?  9  , 

*  c  -  -fr  -  -fr  if  ’  >  <34> 

Equation  34  does  not  adequately  describe  the  motion  of  electrons 
because  it  does  not  include  electron  loss  due  to  pitch  angle  scattering. 

In  most  simple  models,  loss  of  plasma-sheet  particles  by  precipitation  or 
charge  exchange  increases  rapidly  as  the  flux  tube  approaches  the  earth 
and  becomes  smaller  in  volume.  Thus  the  loss  process  tends  to  erode  the 
inner  edge  of  the  plasma  sheet.  We  estimate  the  erosion  of  the  electron 
plasma  sheet  by  keeping  track  of  two  electron  inner  edges  for  each  X,  one 
representing  25%  loss  and  the  other  75%  loss. 

For  strong  pitch  angle  scattering  one  can  estimate  the  precipitation 
life  time  x  as 


7  '  [F'ULF°^iC;ont°ntJ°f  f'UX  *  <W/?W,rl)-l  (35, 

x  (/ds/B)-4/3 

Given  x  we  approximate  the  electron  inner  edge  erosion  by 
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The  inner-edge  motion  for  electrons  is  obtained  by  adding  the  rates  given 
in  (34)  and  (36). 

We  must  now  convert  (34)  and  (36)  to  a  form  suitable  for  numerical  solution, 
and  first  discuss  the  numerical  method  for  stepping  along  in  time.  The 
method  we  have  chosen  is  the  simplest  first-order  one: 

8a:J 

a-i(t+At)  =  aJ(t)  +  — r£— (t)  At  (37) 

k  k  dt 

More  accurate  and  intricate  numerical  methods  for  stepping  along  in  time, 
such  as  second-order  Runge-Kutta  and  a  second-order  predictor-corrector 
method,  were  considered.  The  simpler  scheme  (37)  was  used  because  it  was 
much  easier  to  incorporate  into  our  complex  program,  and  because  analysis 
showed  it  equal  or  superior  to  the  more  complicated  method  with  regard  to 
stability  of  short-wavelength  ripples  of  the  inner  edges.  Flute-type 
ripples  physically  should  decay  exponentially,  especially  rapidly  in 
low-conductivity  regions.  However,  if  the  time  step  is  long  compared  to 
the  physical  decay  time,  the  system  becomes  numerically  unstable,  for  any 
of  the  numerical  methods  we  tested,  and  it  is  the  requirement  of 
stability  against  boundary  ripples  that  limits  our  time  step.  The  more 
complicated  numerical  methods  could  not  improve  on  (37)  with  regard  to 
maximum  allowed  time  step.  For  time  steps  short  enough  to  satisfy  the 
stability  requirement,  the  additional  accuracy  of  the  higher-order 
methods  proved  unimportant. 

We  should  mention  that,  in  order  to  eliminate  a  numerical 

instability  that  tends  to  ripple  inner  edges  on  the  night  side,  we  do  not 

use  (37)  to  compute  aJ(t  +  At)  for  every  j-value.  Instead,  for  even 

k 
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j-values  on  the  night  side,  we  compute  a^(t  +  At)  by  linearly 
interpolating  between  the  a£(t  +  At)  values  for  odd  j. 
p.  Accuracy  tests 

The  complexity  and  the  size  of  our  numerical  system  make  it 
difficult  to  evaluate  the  (error  in  many  conventional  ways.  We  there¬ 
fore  tried  to  estimate  and  minimize  the  error  by  examining  the  dependence 
of  the  results  on  the  values  of  various  numerical-analysis  parameters. 

The  essential  numerical -analysis  parameters  in  our  program  are  the 
following:  A  t,  AS.A^,!^^  ,I-|ayer,  e  and  W.  We  test  accuracy  by  varying 
each  parameter  separately.  The  difference  equations  represent  the 
differential  equations  exactly  in  the  following  limits:  At,Ae,Ai^,e  -*•  0 

and  Knax  "• 

As  mentioned  in  section  IIJC  we  have  used  At  =  30  sec  as  a  constant 
time  step  throughout  the  event.  We  found  that  taking  a  significantly 
longer  time  step  results  in  oscillatory  instability.  Tests  run  earlier 
have  shown  that  choosing  substantially  smaller  At  causes  negligible 
change  in  the  results. 

We  also  have  built-in  means  for  testing  accuracy  by  doubling  A£ 
and  Aij; ,  as  defined  in  equation  11.  We  found  that  doubling  the  grid 
spacing  did  not  change  the  results  in  a  major  way,  although  details  were, 
of  course,  lost. 

The  parameter  k,^,  the  number  of  particle  species  whose  motions  are 
computed  in  detail,  was  not  varied  from  its  chosen  value  21,  because  that 
number  of  species  seemed  to  provide  the  desired  energy  resolution. 
However,  the  program  was  run  for  various  values  of  1 1 ayer *  ’*e*»  various 
numbers  of  subspecies  interpolated  in  between  the  kmax  species  whose 
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motion  we  follow  in  detail  (section  V.C).  Results  for  I,  =  1  and 

layer 

It,...,..*  9  differed  negligibly,  and  we  therefore  chose  I,  =1. 
layer  layer 

The  convergence  criterion  E  defined  by  (28),  was  set  equal  to  103 
volt2.  Test  runs  with  smaller  e -values  produced  insignificantly 
different  results.  The  weighting  factor  W,  defined  by  (26),  was  set 
equal  to  0.5,  close  to  the  largest  value  for  which  the  interation  scheme 
consistently  converged. 

E.  Numerical  problem  near  local  noon 

The  present  numerical  method  has  a  deficiency  near  local  noon. 

Since  we  assign,  for  a  given  particle  species  and  a  given  local  time, 
only  one  i-value  for  the  inner  edge,  we  cannot  have  a  "ring 
current  wrap-up"  which  will  occur  in  reality  when  a  species  of  particles 
attempts  to  form  a  complete  ring  current.  The  problem  becomes  apparent 
near  local  noon  (see  e.g.  Figure  4  in  paper  2)  for  the  high  energy  ions, 
which  are  prevented  from  forming  a  ring.  In  general  this  problem  could 
not  be  overcome  using  the  present  numerical  method  and  we  therefore 
stopped  the  simulation  at  1300  UT.  However  we  did  "prevent"  the 
particles  from  going  to  infinity  near  local  noon  by  holding  them  at  ^ 
10.5  if  the  electric  field  indicated  a  sunward  motion  in  this  region. 
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VI.  OVERVIEW  OF  RESULTS 


A.  Survey 

This  section  is  devoted  to  a  survey  of  numerical  results  from  our 
substorm  simulation.  More  results,  analysis  and  comparison  with  data  are 
presented  in  paper  2  and  paper  3.  (See  also  Harel  et  al. ,  [1979].) 

Four  computer  runs  were  performed  covering  the  four  hour  span  from 
0900  UT  (an  hour  before  onset)  to  1300  UT  (well  into  the  recovery 
phase).  Table  2  summarizes  those  runs. 

The  first  run  has  a  peak  polar-boundary  potential  drop  of  80  kV 
(following  the  solid  line  in  Figure  3).  The  conductivity  model  for  this 
run,  although  including  a  smoothing  function  at  the  plasma's  inner  edge, 
was  less  smooth  than  the  other  three.  Induction 
electric  fields  were  estimated  by  means  of  our  fictitious  substorm 
current  loop  (Figure  7)  turned  on  at  onset  (see  Section  IV. B). 

Run  2  was  similar  to  the  first  run  with  the  exception  of  a  minor 
numerical  method  correction  and  a  greater  smoothing  of  the  conductivity. 
By  comparing  these  results  with  run  1  we  can  learn  the  effect  of  the 
conductivity  smoothing  on  the  numerical  results. 

The  third  run  is  similar  to  the  second  except  that  it  omits  the 
substorm  current  loop.  For  this  run,  the  magnetic  field  model  was  time- 
independent,  in  order  to  evaluate  the  effect  of  the  induction  electric 
field  on  plasma  motion. 

Run  4  was  performed  with  an  estimated  potential  drop  of  140  kV  at 
1150  UT  (dashed  curve  in  Figure  3).  Further  analysis  indicates  that  we 
considerably  overestimated  this  peak  potential  drop. 
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We  believe  run  1  to  be  the  most  realistic  because  of  its  sharper 
conductivity  gradients  and  80  kV  peak  potential  drop;  therefore  most  of 
the  results  presented  here  and  in  the  following  two  papers  were  obtained 
from  that  run. 

In  order  to  be  as  independent  as  possible  of  initial  conditions,  we 
actually  started  all  4  runs  at  0700  UT  (two  hours  prior  to  the  beginning 
of  the  growth  phase).  From  past  experience  we  found  that  in  two  hours, 
magnetospheric  time  and  for  a  50  kV  cross-polar-boundary  potential  drop, 
the  plasma  configuration  and  electric  fields  reaches  near-equilibrium 
almost  regardless  of  initial  conditions  [Hare!  and  Wolf,  1976],  provided 
the  particles  were  initially  put  at  L  £  8.  However,  with  the  present  20 
kV  potential  drop,  the  process  is  slower  and  the  system  progresses  only 
part-way  toward  equilibrium. 

Figure  11  describes  our  assumed  initial  configuration  in  the 
equatorial  plane  of  the  magnetosphere.  Inner  edges  of  3  sample  particle 
species  are  also  plotted  (symmetric  and  initially  placed  at  L  £  10). 
Equipotenti al  lines  6  kV  apart  are  plotted  in  the  rest  frame  of  the 
rotating  earth.  The  electric-field  pattern  is  distorted  because  our 
arbitrarily  assumed  initial  particle  distribution  implied  intense 
localized  partial  ring  currents  and  thus  large  ionospheric  currents  and 
electric  fields. 

Figure  12  shows  the  potential  pattern  at  0900  after  the  initial  2 
hours'  "quiet  time"  polar  boundary  potential  drop.  As  expected,  the 
plasma  sheet,  drifted  in  sunward,  as  described  by  the  3  sample  "inner 
edges"  at  L  x  8.  The  outer  circle  represents  the  mapping  of  our  polar 
boundary  to  the  equatorial  plane.  Note  the  extensive  electric  field 
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shielding  equatorward  of  the  inner  edges,  especially  on  the  nightside; 

(see  also  Jaggi  and  Wolf  [1973]). 

Figure  13  shows  the  electrostatic  potential  distribution  for  0900  as  viewed 
from  above  the  north-polar  ionosphere.  The  equipotential  curves  shown  in 
Figure  11  are  here  plotted  in  the  ionosphere.  Curve  I,  representing  our 
assumed  "polar  boundary" ,  is  an  offset  circle  centered  2  degrees  from  the 
north  pole  toward  midnight.  Shielding  is  obvious  again  equatorward  of 
the  particles'  inner  edges  (  ^  6(f).  The  polar  cap  region  inside  of 
curve  I  (mapped  to  the  region  outside  of  this  curve  in  the  equatorial 
plane  in  Figure  12  and  the  following  figures)  is  not  included  in  our 
model.  Therefore,  electric  field  and  particle  results  shown  in  this 
region  are  merely  extrapolations. 

Figure  14  shows  the  potential  contours  including  the  Earth's  corotation 
potential  for  the  same  universal  time  (0900).  The  earth's  corotation 
field  dominates  for  L  <  5,  represented  by  nearly  circular  equipotentials, 
but  falls  off  at  larger  L  values  where  the  potential  pattern  is  similar 
to  Figure  12.  The  boundary  between  the  two  regions  is  the  "saddle  point 


equipotential"  (termed 


later  an  Alfveil  layer)  of  -15.4  fcV. 


The  next  group  of  figures  is  a  survey  of  some  of  the  results  from 
our  main  run  (run  1  in  table  2).  Figures  15-19  show  equatorial 
potential  contours  in  the  rest  frame  of  the  rotating  earth  and  inner  edge 
locations  at  1000-  eUT,  1000+  e (onset),  1050  (highest  observed 
potential  drop)  and  1300  (partial  recovery)  (Figure  3).  Two  features 


ar  e  iiiusl  obv  iuuS. 


(i)  Injection.  As  time  progresses,  plasma  sheet  particles  are  injected 
sunward  and  around  the  flanks  to  form  an  asymmetric  partial  ring  current 
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at  L  3:5-8  R £  .  Although  others  [e.g. ,  Mcl lwain  [1974],  Roederer  and 
Hones  [1974],  Smith  et  al.  [1979]]  have  reproduced  ring  formations,  their 
models,  unlike  ours,  do  not  include  self-consistent  electric  fields  that 
both  affect  and  are  affected  by  plasma  motion.  As  time  progresses,  the 

inner  edge  of  the  fiigh  energy  ions  drifts  inward  toward  earth  at  dusk 

! 

faster  than  low  energy  ions  and  all  electrons,  giving  rise  to  dispersion 
relations  at  synchronous  orbit  similar  to  those  reported  by  DeForest  and 
Mcl lwain  [1971].  (For  a  more  detailed  discussion,  see  paper  2).  These 
results  point  in  favor  of  the  conventionally  assumed  mechanism  for  the 
formation  of  the  main-phase  ring  current:  the  inner  edge  of  the  plasma 
sheet  is  driven  inward  from  the  tail  to  produce  a  ring  current  with  the 
observed  particle  dispersion.  For  more  discussion  of  this  point,  see 
Wolf  and  Hare!  [1979]. 

( i i )  Shielding.  Although  the  cross-polar-boundary  potential  drop 
increases  to  a  value  of  about  80  kV,  most  of  the  electric  field  is 
concentrated  outside  of  the  most  equatorward  inner  edge  with  some  leakage 
to  the  lower  t  shells.  By  1300  (Figure  19),  shielding  of  the  near-Earth 
region  is  quite  effective.  In  our  case  the  shielding  is  done  mostly  by 
plasma-sheet  ions,  since  most  of  the  electrons  are  lost  by  strong 
pitch-angle  scattering  [Southwood  and  Wolf,  1978],  There  also  seems  to 
be  a  dawn-dusk  asymmetry  in  the  electric  field.  Whereas  on  the  dawn 
side,  electric  fields  decline  smoothly  with  decreasing  latitude,  on  the 
dusk  side,  the  strongest  poleward  electric  fields  generally  occur  well 
equatorward  of  the  polar-cap  boundary.  (For  further  discussion  and 
comparison  of  predicted  and  measured  electric  fields,  see  paper  2). 

Figure  20  shows  the  total  electric  field  in  the  equatorial  plane  at 
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1300.  While  the  corotation  field  dominates  at  L<  4  because  of  the 
earth's  rotation  and  shielding  of  the  potential  electric  field  by  the 
plasma  sheet's  inner  edge,  the  dawn-dusk  convection  field  dominates  at 
large  L  values.  Induction  fields,  although  included,  are  very  small  with 
the  exception  of  the  westward  component  in  the  region  near  local  midnight. 

The  resulting  E  x  B  velocity  flow  field  is  shown  in  Figures  21  and 
22.  The  arrows  show  the  instantaneous  velocity  of  cold  plasma  particles 
in  the  equatorial  plane  (Figure  21)  and  mapped  into  the  ionosphere 
(Figure  22).  At  1800  local  time,  we  can  see  the  flow  reversal  from 
corotation  to  sunward  convection  occurring  at  about  65°  invariant 
lat itude. 

Figures  23-26  show  contours  of  constant  equivalent  potential 

V-V  .  +(x/q) (/ds/B)~2^3  for  various  x’s  and  various  times, 

corotation'  ' 

The  A's  are  "energy  invariants"  of  the  motion  of  various  type  particles 
and  are  defined  in  equation  1.  The  particle  thermal  energy  E  is  related  to 
X  by  equation  1,  which,  for  a  dipole,  can  be  approximated  as 

E  £  (987  eV)U*8/3  (38) 

where  A  has  units  of  eV  r^/ 3^-2/ 3 

In  steady  state,  the  lighter  curves  would  describe  trajectories  of 
the  ions  as  they  drift  around  the  Earth.  One  can  estimate  the  closeness 
of  the  particular  ion  species  to  equilibrium  by  evaluating  the  distance 
from  the  Alfv?n  layer  to  the  ion's  inner  edge  (both  plotted  in  Figures 
23-2G).  Using  this  criterion,  the  higher  energy  ions  (Figures  25-26)  are 
fairly  close  to  "equilibrium"  by  1300  UT,  while  the  lower-energy  ions. 


having  much  slower  drift  speeds,  have  not  come  to  equilibrium  yet 
(Figures  23-24).  For  more  detailed  discussion  see  Paper  2,  Section 
II. B.  The  high  energy  ions  would  have  formed  a  complete  or  nearly 
complete  ring,  corresponding  to  a  recovery-time  ring  current.  They  did 
not  do  this,  however, because  our  numerical  procedure  does  not  allow  an 
inner  edge  to  cross  local  noon  twice  (section  V.E). 

B.  General  Comments 

The  results  shown  in  this  and  the  following  two  papers  represent  a 
"first  try"  at  a  substorm  simulation.  This  is  also  our  first  real 
confrontation  with  data,  and,  as  a  result,  we  have  gained  useful 
experience  from  the  comparison.  We  thus  can  point  to  some  improvements 
in  our  model  inputs  necessary  to  achieve  better  agreement  with  data,  such 
as  better  time  dependent  conductivity  and  magnetic  field  models,  improved 
time-dependent  boundary  potentials,  inclusion  of  a  pre-existing  ring 
current,  a  more  realistic  plasma-sheet  energy  spectrum,  and  an  improved 
formula  for  electron  precipitation.  Abundant  data  are  also  necessary  in 
future  simulations  in  order  to  improve  our  input  parameters. 

General  improvements  in  the  model  should  include  better  numerical 
treatment  of  the  region  near  local  noon,  the  effects  of  neutral  winds, 
field-aligned  potential  drops,  and  allowing  the  system  to  come  to  a  more 
relaxed  state  before  the  beginning  of  the  growth  phase. 

We  believe  that  despite  its  limitations,  our  model  gives  a  great 
deal  of  insight  into  magnetospher ic  dynamics,  and  we  are  encouraged  by 
the  level  of  agreement  between  model  predictions  and  observations. 

The  following  paper  (paper  2)  presents  detailed  comparisons  with 
data  for  the  region  L  ^  6.  Paper  3  deals  with  the  electric  fields  at 
lower  L-values  and  the  resulting  motion  of  the  plasmapause. 
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APPENDIX 


One  purpose  of  this  appendix  is  to  derive  the  formula 

x  §e  x  Ve  [(/ds/B)"2/3] 

-Gee  "  q”  — 


(A.l) 


for  the  average  gradient/curvature  drift  velocity  of  the  equatorial 
crossing  point  of  a  particle  of  charge  q,  for  the  case  of  an  isotropic 
distribution  of  particles;  here  X  is  an  energy  invariant  defined  by 
equation  1  (section  II .A),  =  equatorial  magnetic  field,  and  = 

2-dimensional  equatorial  gradient  operator.  This  formula  has  been 
derived  by  others  [e.g.,  Vasyl iunas,  1968;  Southwood,  private 
communication],  but  no  derivation  has  been  published,  as  far  as  we  know. 
The  derivation  presented  here  arose  out  of  conversations  with  D.  J. 
Southwood  and  G.  M.  Erickson. 

The  second  purpose  of  this  appendix  is  to  show  that,  for  particles 
drifting  according  to  a  law  of  the  form  (A.l),  the  energy  invariant  X, 
and  the  number  of  particles  n  per  unit  magnetic  flux  are  constant  along  a 
drift  path. 

Consider  a  symmetric  magnetospheric  magnetic-field  configuration, 
with  an  equatorial  plane  that  is  always  normal  to  the  local  northward- 
directed  magnetic  field.  Consider  only  field  lines  that  go  through  the 
earth,  where  the  magnetic  field  is  very  strong.  Assume  that  the  value  of 
the  particle  kinetic  energy  E  (gyration  and  bounce  motion)  is  uniquely 
determined  by  the  equatorial  crossing  point  of  the  orbit  xe,  and  by  some 
other  parameters  that  are  constant  along  a  drift  path.  (In  the  case  of 


54 


adiabatic  motion  these  other  parameters  might  be  the  particle's  first  two 
adiabatic  invariants.)  In  the  case  where  a  pitch-angle  scattering 
process  keeps  the  distribution  function  isotropic,  but  the  scattering 
process  involves  no  appreciable  energy  transfer,  the  appropriate 
invariant  is  A,  where  , 

E  (A,xe)  =  A  [H(xe)]"2/3  (A.2) 

where  E  (xg)  represents  an  effective  volume  occupied  by  the  particles. 

We  will  derive  an  explicit  expression  for  E  (xg)  later  (equation  A. 10). 
Equation  (A. 2)  is  just  the  statement  that,  in  adiabatic  compression  of  an 
ideal  monatomic  gas,  the  mean  thermal  energy  is  proportional  to  the  2/3 
power  of  the  density. 

If  the  particles  are  E  x  B  and  gradient/curvature  drifting  in  a 
static  magnetic-field  configuration  of  the  type  described  above,  and 
drift  velocities  are  small  compared  to  thermal,  then  the  conservation- 
of-energy  relation  can  be  written 

(*EBe  *  ~Gce)  *Ve  EqV{?e)  +  E(X’*e)]  =  0  {A*3) 

where 

~EBe  E  r?  <-V>  x  5e  <A-4> 

Be 

is  the  E  x  B-drift  velocity  in  the  equatorial  plane,  and  V  =  electro¬ 
static  potential.  Using  the  fact  that  v^  cannot  depend  on  V,  the  fact 
that  (A. 3)  holds  for  in  an  arbitrary  direction  in  the  equatorial 
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plane,  and  some  vector  identities,  we  obtain  from  (A. 3): 


-Gee  '-<l8e2rl  V<A.Xe)  *  ?, 


(A.5) 


Any  collection  of  particles  that  drifts  according  to  a  law  of  the 


Eq  x  B  V  Y  x  B 

v  =  ~e  +  _§ _ I®. 

~e  2  2 

Bo  B  * 

e  e 


(A. 6) 


conserves  n,  the  number  of  particles  p°r  unit  magnetic  flux.  To  show 
this,  we  write  the  law  of  conservation  of  particles,  mapped  to  the 
equatorial  plane,  in  the  form 


3(nB  ) 

Hstr5- +  * 0 


(A. 7) 


where  nBg  =  particles/area  mapped  to  the  equatorial  plane,  and  we  have 
assumed  that  no  particles  are  lost.  We  now  let  z  =  northward  unit 
vector  in  equatorial  plane  and  use  (A. 6)  to  rewrite  (A. 7)  in  the  form 


Be(lt  +  VVe>  *  +  1  <  jf  +  V[(§e  +  V°  x  ^  >  =  0 


(A. 8) 


or,  using  vector  identities  and  Faraday's  law  (7  x  E  =  -38/3t), 


<5t  *  Vve>  "  =  0 


(A. 9) 


We  have  shown  that  the  convective  derivative  of  n  is  equal  to  zero,  which 
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was  one  of  the  purposes  of  this  appendix. 

Equation  A. 9  tells  us  how  the  effective  volume  5  (xg)  changes  in 
this  convection.  Namely,  if  a  set  of  particles  initially  occupied  one 
unit  of  magnetic  flux,  they  will  forever  occupy  one  unit  of  flux.  The 
volume  corresponding  to  one  unit  of  magnetic  flux  is  thus  / ds/B,  where 
the  integral  covers  the  entire  flux  tube,  and  this  is  the  effective 
volume  in  (A. 2).  Thus  the  relationship  between  the  kinetic  energy  E  of 
an  individual  particle  and  the  energy  invariant  X  is 

E  =  X  {/ds/B]"Z/3  (A.  10) 

Here  and  elsewhere  in  the  paper,  the  limits  of  /ds/B  are  assumed  to  be  in 
the  northern  and  southern  ionospheres.  Substituting  (A. 10)  in  (A. 5) 
yields  (A.l),  which  is  the  equation  we  set  out  to  derive. 
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Table  1.  Model  assumptions. 


A.  Characteristics  of  the  model 

(1)  Region  of  model  =  inner  magnetosphere  and  ionosphere;  l£  10, 

invariant  latitute  -  70  .  ■ 

(2)  Olson-Pf itzer  analytic  magnetic  field  model  and  substorm  current 
loop  (not  self-consistent). 

(3)  Time  dependent  ionospheric  conductivity  model,  including 
day-night  asymmetry  and  auroral  enhancement. 

(4)  V  •  j  =0  both  in  magnetosphere  and  ionosphere. 

n,' 

Self-consistently  calculated  current  system,  consisting  of 
horizontal  ionospheric  currents,  magnetospheric  ring  currents, 
and  Birkeland  currents  connecting  the  two. 

(5)  E  =  _  v  Vin  ionosphere. 

(6)  Isotropic  distribution  of  particle  pitch  angles  in  the 
magnetosphere. 

(7)  Discrete  particle  energy  spectrum.  Electron  energy  ^1-4  keV; 
ion  energy  >\,  500  ev  -  60  keV. 

(8)  Electron  loss  via  strong  pitch-angle  scattering. 

(9)  Gradient,  curvature  and  E  x  B  drifts  included.  Polarization 

Or  'V* 

currents  excluded. 

( 10) Time-independent  particle  input  at  the  high-L  boundary  of  the 
model.  Maxwellian  input  distributions  are  assumed  with  Tj  % 

4.5  keV,  Tq  £  1.5  keV  at  L  X  10. 

B.  Not  Included  in  the  Model 

(1)  Field-aligned  electric  fields. 
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(2)  Ionospheric  neutral  winds. 

(3)  Pre-existing  ring  current.  All  particles  assumed  to  originate 
at  L  £  8. 

(4)  Polar-cap  phenomena.  Solar  wind  and  polar  cap  phenomena  enter 
model  only  through  boundary  conditions. 

(5)  Equatorial  electrojet.  Equatorial  boundary  condition  is 


latitude. 
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Table  2.  Computer-Simulation  Runs 


Computer  Run 

Peak  Polar-Boundary 
Potential  Orop 

Conductivity 

Model 

Induction 
Electric  Field 

1 

80  kV 

minimum  smoothing 

yes 

2 

80  kV 

greater  smoothing 

yes 

3 

80  kV 

greater  smoothing 

no 

4 

140  kV 

greater  smoothing 

yes 
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Figure  Captions 

Figure  1.  Corresponding  regions  in  the  ionosphere  and  the  equatorial 
plane;  mapping  is  done  via  the  pre-substorm  magnetic-field 
model.  The  circle  I  in  the  ionosphere  represents  the  assumed 
polar  boundary  and  maps  to  curve  I  in  the  equatorial  plane. 

Our  computer  model  applies  to  the  region  equatorward  and 
earthward  of  curve  I. 

Figure  2.  Overall  logic  diagram  of  our  program.  Arrows  indicate  flow  of 
information  in  the  program.  Dashed  lines  indicate  features 
that  we  plan  to  incorporate  in  the  program,  but  are  not 
included  yet.  The  rectangles  at  the  corners  of  the  central 
pentagon  represent  basic  parameters  computed.  Input  models 
are  indicated  by  rectangles  with  round  corners;  input  data  are 
indicated  by  curly  brackets.  The  program  cycles  through  the 
entire  main  loop  (including  all  the  rectangular  boxes). every 
time  step  At  (approximately  every  30  seconds).  The  basic 
equations  that  the  computer  uses  or  solves  are  described 
briefly  by  words  or  symbols  next  to  the  logic-flow  lines. 

Figure  3.  Fort  Churchill  H-magnetogram  and  polar  cap  potential  drop  for 
19  September  1976.  Top  panel  shows  the  negative  bay  for  the 
event  we  model  (onset  at  about  1000  GMT,  maximum  at  1050). 

The  1050  dip  is  £520y  .  A  small  disturbance  is  seen  around 
0600  GMT.  Bottom  panel  gives  estimates  of  polar  cap  potential 
drop  (  /£.  cH)  across  the  polar  cap  for  various  crossings  by 
satellite  S3-?.  Rectangular  boxes  represent  potential  drops 
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from  field  reversal  to  field  reversal,  the  sizes  of  which 
gives  the  uncertainties  in  the  integration.  Vertical  lines 
(between  1000  UT  and  1200  UT)  represent  our  polar  boundary 
potential  drop,  where  we  use  the  region-1  currents  to  locate 
the  boundary.  The  solid  line  represents  our  main  choice  runs 
(  1,  2,  and  3  in  Table  2.)  Dashed  curves  give  an 

alternative  choice  with  peak  potential  drop  of  140  kV  (run  4). 

Figure  4.  Precipitating  electron  energy  flux,  east -west  magnetic 

perturbations,  a  B,  and  forward  component  of  electric  field 
for  orbit  4079A- Southern  hemisphere.  (S3-2  satellite  data). 
Arrows  in  bottom  panel  point  to  the  most  equatorward  electric 
field  reversal;  vertical  lines  show  the  equatorward  edge  of 
the  region-1  currents. 

Figure  5.  Same  as  figure  4  except  that  data  from  orbit  4079A  North  is 
plotted  (polar  cap  crossing  of  roughly  1050  UT).  Also  shown 
on  Figure  5  are  features  of  the  Lockheed  particle  detector  on 
the  S3-3  satellite  at  comparable  invariant  latitudes  (nearly 
simultaneous):  the  top  panel  shows  features  of  the  electron 
detector  (weak,  hot,  cool),  plus  the  locations  of  "inverted 
V's"  (shown  by  carats).  The  middle  panel  shows  features  of 
the  ion  fluxes  (weak,  strong),  plus  the  locations  of  upward 
flowing  ions  (indicated  by  *’  UF  I  ” ) . 

Figure  6.  Same  as  Figure  4.  South  pole  crossing  occurring  at  about  1140 
UT. 

Figure  7.  Schematic  diagram  of  the  substorm  current  loop.  View  is  from 
above  the  equatorial  plane.  Arrows  indicate  direction  of  the 


current.  Loop  includes  an  equatorial  segment  (eastward 
current  in  the  tail),  a  westward  electrojet  in  the  polar 
ionosphere  and  two  segments  of  field  aligned  currents. 

Figure  8.  Contours  of  constant  magnetic  field  in  the  equatorial  plane. 

View  is  from  above  the  north  magnetic  pole.'  Spacing  between 
successive  contours  follows  a  logarithmic  scale.  Values  (in 
Gammas)  are  given  for  a  few  contours. 

Figure  9.  Pedersen  and  Hall  conductivity  profiles  for  various  latitudes 
in  the  ionosphere.  Height  integrated  conductivity  in  mhos  is 
plotted  vs.  local  time  for  constant  invariant  latitudes  of  55° 

,  6CP  ,  65°  and  7(F  ;  and  constant  universal  time  of  1150 
(peak  polar-cap  potential  drop).  Tic  marks  are  spaced  2  mhos 
apart.  The  values  correspond  to  run  1  ,  which  features  the 
least  smoothing  at  the  plasma's  low-latitude  edge  of  the 
auroral  zone. 

Figure  10.  Auroral  conductivity  profiles  for  various  times  during  the 
substorm.  Plot  format  is  the  same  as  Figure  9.  Latitude  is 
the  same  for  all  profiles  (7Cf  =  poleward  region  of  the 
auroral  oval).  Conductivity  profiles  are  given,  starting  from 
the  top  panel,  for  0900  (beginning  of  growth  phase),  1000  +  e 
(onset),  1050  (peak  of  substorm),  1150  and  1300  (at  which  we 
terminated  our  runs). 

Figure  11.  Initial  potential  distribution  and  plasma  "inner  edge" 

locations,  in  the  equatorial  plane.  The  sun  is  to  the  left. 
Fquipotentials  are  6  kV  apart  and  plotted  in  the  rest  frame  of 
the  Farth.  Also  plotted  are  "inner  edges"  of  three  sample 
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particles  (out  of  21  different  "energy  invariant"  species) 
that  we  monitor.  The  inner  edge  was  arbitrarily  chosen  to 
be  £  10  for  this  initial  configuration  (corresponds  to  0700 
UT). 

Figure  12.  Potential  contours  without  corotation  in  the  equatorial 

plane.  Potential  distribution  and  locations  of  the  "inner 
edges"  were  computed  two  hours  after  the  initial  arbitrary 
starting  condition  shown  in  Fig.  11.  The  closed  solid  curve 
at  L^IO  is  the  boundary  of  our  calculation. 

Figure  13.  Equipotenti al  contours  in  the  ionosphere.  The  view  is  from 
above  the  north  pole.  The  spacing  between  two  successive 
contours  is  6  kV.  The  inner  circle  is  our  assumed  boundary 
curve  I  and  is  offset  2 °  toward  midnight.  Also  plotted  are 
the  inner  edges  of  3  of  the  21  types  of  particles  we  model. 

Figure  14.  Same  as  Figure  12  except  for  the  inclusion  of  the  corotation 
potential.  The  saddle  point  potential  (15.4  kV  contour)  is 
also  plotted. 

Figure  15.  Same  as  Figure  12  except  for  a  different  time  1000-e  UT  (just 
before  onset). 

Figure  16.  Same  as  Figures  12  and  15.  Equipotentials  are  plotted  for 
1000+e  UT,  following  the  conductivity  jump  at  onset. 

Figure  17.  Same  as  Figures  12-16.  Contours  are  plotted  for  1050  UT  (peak 
of  expansion  phase). 

Figure  18.  Contours  of  electrostatic  potentials  without  corotation  in  the 
equatorial  plane  at  1150  UT  (corresponding  to  peak  in 
cross-polar-cap  potential  drop).  Format  is  the  same  as 


previous  figures. 

Figure  19.  Same  as  Figure  17,  but  for  1300  UT,  deep  into  the  recovery 

phase.  All  of  our  simulation  runs  were  terminated  at  1300  UT. 
Figure  20.  Total  electric  field  in  the  equatorial  plane  (1300  UT). 

Arrows  give  the  direction  and  magnitude  of  the  combined 
potential  and  induction  field. 

Figure  21.  Velocity  flow  field  in  the  equatorial  plane.  This  velocity 

field  is  E  x  B/82drift  velocity  where  E  is  the  electric  field 
^  %  % 

of  Figure  19.  Arrows  give  the  instantaneous  velocity  vector 
of  cold  plasma  particles. 

Figure  22.  Total  flow  field  in  the  ionosphere  for  1300  UT.  Noon  is 

upward.  Circles  representing  constant  latitude  are  spaced  10° 
apart. 

Figure  23.  Instantaneous  trajectories  and  inner  edges  (dashed curves) 

forionswith  A  =  478  eV  R^^y-2^3  .  The  Alfven-layer  and  the 
inner  edge  curves,  for  a  given  time,  are  darker  than  the  other 
trajectories.  The  plots  are  equatorial -plane  views  for  various 
phases  of  the  substorm. 

Figure  24.  Same  as  Figure  23,  but  for  a  =  1730  eV  R^\”^3 

Figure  25.  Same  as  Figures  23-24,  but  forA  =  3880  eV  R£2/^y~2^3 

Figure  26.  Same  as  Figures  23-25,  but  for  y  =  8650  eV  R^3  Y 
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